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The kinematic properties of unsteady highly non-linear 3D wave groups have been inves¬ 
tigated using a numerical wave tank. Although carrier wave speeds based on zero-crossing 
analysis remain within ±7% of linear theory predictions, crests and troughs locally un¬ 
dertake a systematic cyclical leaning from forward to backward as the crests/troughs 
transition through their maximum amplitude. Consequently, both crests and troughs 
slow down by approximately 15% of the linear velocity, in sharp contrast to the predic¬ 
tions of finite amplitude Stokes steady wavetrain theory. Velocity proHles under the crest 
maximum have been investigated and surface values in excess of 1.8 times the equivalent 
Stokes velocity can be observed. Equipartitioning between depth-integrated kinetic and 
potential energy holds globally on the scale of the wave group. However, equipartitioning 
does not occur at crests and troughs (even for low amplitude Stokes waves), where the 
local ratio of potential to total energy varies systemically as a function of wave steepness 
about a mean value of 0.67. 


1. Introduction 


Wind-generated ocean waves at the air-sea interface propagate predominantly in group 
or packet structures. Group structure behaviour is not limited to ocean waves, but occurs 
naturally in many systems. These unsteady wave groups can exhibit a complex nonlinear 
life cycle, especially in focal zones where there is rapid wave energy concentration. Their 
structure and the associated carrier wave propagation properties of dispersion, direction¬ 
ality and nonlinearity, and their interplay, represents a knowledge gap beyond present 
analytical treatment. Understanding natural ocean wave propagation is not only relevant 
academically but is potentially important when studying ocean-atmosphere exchanges, 
as well as the design of structures such as open ocean platforms. Crest speeds are a 
key consideration in the assessment of wave impact loading on structures, where flow- 
induced forces are proportional to the square of the water velocity. In other applications, 
measurements of in itial whi t ecap speeds have been proposed as a method for inferring 
ocean wavelengths ( PhillipsI ( 1985 1) assuming that the Stokes dispersion relationship is 
applicable. Any systematic crest speed slowdown just prior to breaking onset will signif¬ 
icantly alias wavelength determinations from the quadratic relationship between Stokes 
wavelength and speed. 

Stokes classical deep water wave theory (( Stoked ( 18471) ) was developed for a steady, 
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uniform train of two-dimensional (2D) non-linear, deep-water waves of small-to-intermediate 
mean steepness ak (= 27r^), where a is the wave amplitude and A is the waveleng t h. Th e 
wave speed c increases with ak as found to 5th order in wave steepness by iFentonI (119851) : 


c = Co ( 1 -I- - {ak) + 


— (ak)'^ -(-higher order terms in {ak)^^ 


( 1 . 1 ) 


where cq is the wave speed of linear (infinitesimally steep) waves. 

Increasing the steepness to the Stokes limit {ak ^ 0.42) in eq (11.11) gives the lim¬ 
iting speed of maximally-steep steady waves of I.Icq. Increased wave steepness has 
long been associated with higher wave speeds. Historicall y, Stokes theory has b een the 
defaul t theoretical approach for des cribing ocean wav es (|Kinsmanl (|l9651) and IWiegel 
1 1964l) l. Contemporary examples are Beva et al. ( 2012ll who found close agreement be¬ 
tween monochromatic wave near-surface velocities and pred ictions based on 5th order 
Stokes theory. The observational study of iGrue et ali (|2003l) using PIV to measure ve¬ 
locities under extreme and breaking crests sh owed that invi s cid pr edictions are in good 
agreement with measured particle velocities. Ijohannesseni ( 201011 showed that at the 
crest, near-surface velocities agree with second-order Stokes theory, which is capable of 
describing the free surface kinematics very accurately provided that the local underlying 
regime of free waves can be identihed. A principal constraint of the Stokes wavetrain 
framework is the imposed uniform spatial and temporal periodicity, which provides the 
basis for the extensive analytic treatment in the literature of this intrinsically nonlinear 
free surface problem. 

However, ocean waves propagate naturally in unsteady groups that evolve dynamically. 
Relaxing the periodicity constraints all ows addi t ional degrees of freedom to develop as 
shown in a suite of laboratory studies. iMelvillel ( 19831) fo und a wave spe ed variation of 
between -17% and -|-32% around the linear phase velocity. Shemer ( 2013h demonstrated 
that for a wave group with a broad spectral bandwidth, the crest propagation velocity 
of the dominant waves differs signihca ntly from both th e ir ph a se and group velocities. 
Deep- w ater breaking wave studie s (e.g. Rapp &: Melville! ( 1990ll. Stansell fc MacFarlanj 
( 2002ll . Ijessup fc Phadni^ ( 20051) . [Melville fc MatusoX ( 2002 1. iMelville et al\ ( 2002ll ) all 
found a significant (0(20%)) breaking crest speed slowdo wn relative to the expected 
line ar phase velocity (Eq. 11.11). For focusing wave packets, Johannessen fc Swan ( 200ll) 
and [johannessen fc Swan! (|20()3l) reported a slowdown in the expected crest velocity at 
focus of around 10%. 

Understanding wave crest slowdown behaviour is central to both refining present 
knowledge on water-wave propagation and dy namics, and eff e ctive implementation o f 
Phillips’ spectral f ramew ork for breaking waves ( Phillips! ( 1985 ), Gemmrich et all ( 2008 ), 
iKleiss fc Melville ( 2010l) i. The groupines s of natural ocean wa ves gives rise to a suite of 
differences from Stokes waves predictions. iBaldock et al\ ( 1996h report an ensemble of di¬ 
rectional wave modes of different frequency components focusing in space and time to pro¬ 
duce an unsteady wave-group in the laboratory. Surface eleva tions and subsurface particle 
kinem atics were compared with linear wave theory and the iLonguet-Higgins fc Stewart 
( 19601) second-order solution of the wave-wave interactions. Their study shows that non¬ 
linear wave-wave interactions produce a highly nonlinear wave-group with crests higher 
and t roughs shallower tha n expected from numerical simulation s by Longuet-Higgin^ 
( 19871) and experiments by Miller et al. { 199ll ) . [Sutherland et al. ( 1995 ) previously ob¬ 
tained very similar c haracterisations a l thoug h they also conclude d that the w ave kinemat¬ 


ics were Stokes-like. Song fc Banner ( 2002I) . Banner fc Peirson ( 2007 ) and [Viotti et al. 


(|2014i) report complementary numerical and experimental studies that demonstrate that 
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intra-group and inter-wave energy transfers cannot be neglected when characterising wave 
group kinematics and dynamics. 

In this contribution we show that Stokes characterisations have significant limitations 
when applied to fully non-linea r unsteady wave groups. Specifically, im portant proper¬ 
ties derived from IStokesI ( 18471) theory and subsequent refinements (e.g. Fenton ( 1985ll l 
change fundamentally when stationarity is relaxed. Under wave crest maxima, the veloc¬ 
ity profiles, energetics and energy partitioning all differ syste matically from those of a 
comparably steep Stokes wave. The unsteady leaning of waves ( TavfunI ( 198^ 1 is shown 


to be a critical aspect, determining the actual wave crest speed. Our preliminary nu¬ 
merical studies of chirped wave packets have alrea dy been verified agai ns t obser v ations 
of wave groups in the laboratory and in the field ( Banner et all ( 20l4) l. iFedelel ( 2014 1 
provided a theoretical explanation and further validation of the crest slowdown in terms 
of geometric phases. 


2. Numerical simulation 

For mod ulated wav e s, the nonlinear Schrodinger (NLS) equation and its higher-order 
extensions ( Zakharo^ ( 1968 11 have been widely applied. NLS models have been used 
extensively to study the evolution of quasi-monochromatic waves and instabil i ties d ue 
to resonant interactions (e.g. the Benjamin-Feir instability ( Beniamin fc Feiil ( 1967ll l. 
Extensions of the NLS equation include the Zakharov equation which has been applied 
to describe the long-time evolution of the spectrum of weakly nonlinear, dispersive waves. 

There has been a growing interest in the development of three-dimensional models 
which inherently incorporate non-linearity and associated dispersion effects. The broad¬ 
bandedness in both frequ ency and direction of r eal sea states poses significant challenges 
in numerical simulation. Bateman et all ( 2001 1 demonstrated the importance of direc¬ 
tionality and the consequent benefits of efficient wave mo delling during compa r ison o f 
numerical simulations with the laboratory observations of l.lohannessen fc SwanI (j200l[l . 
High-order spectral expansion appr oaches using effi c ient FFT solvers for app lication to 
3D waves have been developed (e.g. Ducrozet et al. ( 2012h. Dias et all ( 201511 1. Another 
option is to solve the full Navier-Stokes equations ( Park et a7r ( 200^ . but viscous flow 
solvers tend to be too dissipative and computationally time-consuming. 

Numerical models of potential 3D wave propagation can be divided into thre e main cat- _ 

egorie s : (al boundary element in t egral methods (BEM): e.g. Bakei_eta^ (Il982h. Bateman_et_aL 
( 2001 1. Clamond fc Grud ( 2001 1. Grilli et al\ (12001 1. Fochesato et al. ( 2007ll . Guvenne fc Grilli 
( 200(3 1. Xue et al. ( 2001 1. Hou fc Zhand ( 200^ Fructus et all 1 20()5ll: (b) hnite element 


method (FEMl e.g. iMa et al. (2001 1; (cl spec t ral methods: e. g. Dommermuth fc Yud 


( 1987l l. West et al. ( 1987 1. Craig fc Siileml ( 1993 1. Niched ( 1998ll . Bateman et al. ( 200 il l. 


Spectral methods based on perturbation expansions are known to be very efficient. 
These methods reduce the water wave problem from one posed inside the entire fluid 
domain to one posed on the boundary alone, thus reducing the dimension of the formula¬ 
tion. This reduction can be accomplished by using integral equations over the boundary of 
the domain (so-called boundary integral methods) or by introducing bou ndary quantities 
which can be expanded as Taylor series for reference domai n geometrie s ( Xu fc Cliivennel 
( 2009l ll. Both approaches have been summarised recently bv iMal ( 20ldl . BEM techniques 
are efficie nt for representing wave prop agation and overturning until the wave surface re¬ 
connects ( Grilli fc: Subramanval ( 19961) 1. 

The present study used a boundary element nu merical wave tank code called WSIM, 
which is a 3D extension of the 2D code developed by Grilli et~^ (198^ to solve the single¬ 
phase wave motions of a perfect fluid. It has been applied extensively to the solution of 


























































































































































4 



finite amplitude wave propagation and wave breaking problems (see chapter 3 of iMa 
( 201011 1. The ideal fluid assumption makes WSIM unable to simulate breaking impact 
subsequent to surface reconnection. However, its potential theory formulation enables 
it to simulate wave propagation in a CPU-efficient way, without the diffusion issues 
of viscous numerical codes. The simulation of wave generation and development of the 
onset of breaking event s can be carried out wi th great precision, as demonstrated by 
Fochesato fc Dias ( 2006ll . Fochesato et al. ( 2007 1. 


vv oiivi nas oeen vaiiaaiea extensively ana snows exce 
(1989ll. Crilli & Svendsen (19901. Grilli & Subramanva 

iieiii energy conservation K mill ei ai. 
(Il9941. Grilli & Subramanval (Il99fill. 

Grilli & Horrillol ( 1997ll. ICrilli et all (l200lh. iFochesato 

(12004). Fochesato & Diad (20061. 


Fochesato et aLl ( 2007ll l. Its kinematical accuracy has been demo nstrated b y com parison 


with the analytical solutions for infinitesimal sine waves found in Phillipd ( 1977 1. 


2.1. Governing equations 

WSIM uses a mixed Eulerian-Lagrangian time-updating scheme for irrotational motion 
described by the velocity potential </> (x, t), in a Cartesian coordinate system x = (x, y, z) 
with constant pressure at the open water surface, z is the vertical upward direction and 
z = 0 the still water surface (Figure [J). 

Fluid velocity is defined as u = Vcj) = {u,v,w). The continuity equation leads to the 
Laplace equation for the potential within the fluid domain U (<) (12.11) . The symbols F and 
F FS are used to denote the entire domain boundary and the free surface respectively. 


= 0, on U 

(2.1) 

Dr 

— = u = V0, on Tps 

(2.2) 

— = -gz p -\/(j)V(j) -, on Fj^s 

Dt 1 p 

(2.3) 

dn(j) = 0, on T\Tfs 

(2.4) 


where r is the position vector of a fluid particle on the free surface, g the gravitational 
acceleration, Pq the atmospheric pressure, p the fluid density and -^ (= ^ + ■ V) 

the Lagrangian (or material) time derivative. 

The free surface (domain Tps) is described by fully-nonlinear kinematic lea 12.21) and 
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dynamic p.3l) equations. Recent developments have implemented a 3D snake wave paddle 
at one vertical face of the domain which is described subsequently. The far face of the 
numerical wave tank from t he paddle has an absorb ing beach which damps any incident 
wave energy as described in Grilli fc Horrillol 1 1997ll . All remaining faces of the domain 
have a zero-flux boundary condition (r\r ps) lea 12.41) . 


2.2. Discretisation 


Boundary geometry and field variables are represented by 16-node quadrilateral elements, 
providing bi-cubic local interpolation of the solution between the nodes. High-order tan¬ 
gential derivatives required for the time updating are calculated in a local curvilinear 
coordinate system, using 25-node fourth-order quadrilateral elements to retain third or¬ 
der accuracy in the spatial derivative. This results in a high-order integration scheme 
over the discretisation e lements. The solver uses a fast PTMA scheme introduced by 
Fochesato k. Diaj ( 2006ll . 


Once preliminary testing of WSIM was complete, two spatial resolutions were used 
during the present simulations to quantify sensitivity to the discretisation. These were 8 
points per wavelength in (so-called) medium resolution and 16 points per wavelength in 
high resolution in the propogation (x) direction. The number of points in the transverse 
y and depth z directions was adjusted to keep element aspect ratios close to 1. 

Even if WSIM internally works with a different reference system, the results in this 
paper have been chosen to be expressed in deep water units. Model time scales are 
non-dimensionalised by the base frequency cUp of the generating paddle signal, and its 
corresponding length scale is the deep water wavelength Ap. These choices impose a 
gravitational acceleration of gnw = 27r and a reference linear velocity Cpw = = 1- 

Non-dimensional lengths can be computed using Lpw = ^ and the time related 
quantities using Tpw = w- 

UJp 

The numerical scheme is a mixed Eulerian-Lagrangian method. Consequently, the mesh 
is subject to Lagrangian drift during the passage of steady, steep wave g roups, especially 
near t he free-surface. For the transitory Class 3 wave groups defined in ISong fc Bannen 
( 2 OO 2 II and used during this present investigation, the mesh did not deform significantly. 


2.3. Wave group generation 

The WSIM wave-maker was chosen as a bottom-hinged flap paddle to generate deep¬ 
water waves. 

The wave-maker motion requires specification of the velocity potential on the boundary 
as follows: 


d(j) 

dn 


= rgU 

bound. 


d'^cj) 

dtdn 


= rgO. + a 

bound. 


d'^4> 


d(j) 


(2.5) 

( 2 . 6 ) 


where Vg is the distance of the centre of the rotation and a is the angle of rotation of the 
wave-maker. 

To ensure smoothness of the boundary conditions and to avoid numerical instabili¬ 
ties at t = 0 created by the moving elements, the paddle oscillation motion 
is damped by a function D (t) which changes smoothl y from 0 to 1, as detailed in 
Crilli fc Subramanval ( 1996l l and ICrilli fc Horrillol (199^. Hence, the boundary condi¬ 


tion can be written as follows: 




































6 


X. Barthelemy, M.L. Banner, W.L. Peirson, F. Dias, M. Allis 

z 


I 



Figure 2. 2D Wavemaker 


d(j) 

dn 

dtdn 


bound. 


bound. 


1 

II 

(2.7) 

= - ^Yx3^D{t) 

(2.8) 


Because of finite tank-length considerations, we chose to investigate the chirped wave 
packet (Class 3) and investigated cases comprising five, seven and nine carrier waves in 
the initial time wave packet. 


2.3.1. Class 3 wave groups 

Using the prescribed ’chirp’ motion from ISong fc Banner ( 2002 ). the following expres¬ 
sion defines Class 3 wave packets, N being the number of waves in the time series, ojp is 
the paddle frequency and 9{t) = k [xui — ct] — 9o is the phase. Ct 2 is the chirp rate of the 
linear modulation: 



1 -I- tanh 


y Ntt 


sin 



/ _ UJpCt2t^ 

V 2 


^1 — tanh 


/ A{(jjpt — 2 Ntt)'\\ 

iv )) 


^ “b {Xconv: 



The paddle motion has an oscillating part moderated by a damping function. To avoid 
discontinuities and assure a smooth numerical start in the code, there is a negative shift 
in time of 5 seconds, to ensure that any initial oscillations are of the order of the double 
precision zero. 


2.3.2. 3D converging wave group 

If 3D waves are to be produced by the snake paddle, the phase $ becomes a function 

‘-coruJi Yconvl 


of {Xrrmjt.Yr.nnv) which sp c cify the lin e ar con vergence point coordinates, as defined by 
Dalrvmple fc Kirb^ ( 198^ : [Palrymold (198^: 
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/I , y ^conv 

Un = arctan- 

T^conv 

^ — k ' y ' sin On ~\- k (^-^conv cos On “t“ Yconv sin On') 

3. Post-processing 

3.1. Determining near-surface guantities: 

A significant challenge within this investigation was determining quantities close to the 
highly-curved free surface. Three complementary solutions have been implemented. Eval¬ 
uating the interior fields is known to be a ’quasi-singular’ problem, because it evaluates 
integration with kernels of the type , r being the distance between the point where the 
field is calculated and a point on the boundary and a is a positive exponent. The closer 
the evaluation points are to the boundary, the larger these terms become and they can 
possibly grow without bound. 


3.2. Interior fields 

The boundary element formulation of 3D potential problems implies: 


a(Xs)(/)(Xs) = /" {q(j>* - dS 
Js 


where Xg is the source point, (fCx-s) is the potential; and g (x) = — is the derivative 

on 

of u along the unit outward normal n at x on the boundary S. S is the boundary of 
the region V of interest and boundary conditions concerning (j) and q are specified on S. 
a (xg) = 1 when Xg G E and a (xg) = i when Xg G S and S is smooth at Xg 
The fundamental solutions (j)* and q* are defined by: 


(x,Xg) =-—, 

47rr 


X,Xs 


(Pn) 

dTTT-S 


where r = x — Xg and r = |r|. 

The flux at a point Xg € E is given by the potential gradient: 


d(j) 

dxs 




dS 


(3.1) 


where 


d(j>* V dq* 1 / n (!■,») 

9xg Awr^ ’ i9xg 47r \ r® 


Equations are discretised on the boundary S by boundary elements Se{e = 1,N) 
defined by the interpolation function. The integral kernels of these equations become 
nearly singular when the distance d — (r, n) between Xg and Se is small compared to the 
size of Se- In the following, we will denote Sg by S for brevity. 

Reconstructing the inner field values relies on the ability to correctly estimate the sum 
iH) over the entire domain. Three complementary methods are used depending on the 
distance d between Xg and Se- 

When the distance d is large enough (relative to the element size), the integration 
on the individual element is carried out by a Riemann quadrature on a classical Gauss- 
Lobatto point distribution. 

When the distance d becomes small, the classical quadrature does not retain enough 
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precision, and the TellesI (1987) method is used. The Telles method consists of binary 
subdivisions of the integration space to maximize precision where the singularity begins 
to manifest itself. The precisi on of this method is a c cepta ble at moderate distance from 
the boundary, as described in iGrilli &: Subramanval (119941) . 

However, the Telles method becomes inefficient when the near-boundary singularities 
are too strong, specifically adjacent to the highly curved surface of steep waves. Conse¬ 
quently, a third method was developed and imp l emen t ed called Proje c tion and Angular _ 

and R a dial Transformat i on (PART) ('see Havamil ( 1990l) , Havami ( 1991 ), Havami fc Matsumotd 
( 1994l) .( Havami ( 2005 q ). Havami ( 200561) ). A change of variables is used to calculate the 
integration on a closed element. 

When applied to the curved quadrilateral element S, for the nearly-singular integral 
over a curved boundary element: 


/ = Us (3.2) 

A sequence of steps is needed in order to effect a change of variables and integrate: 

• 1. Source projection: Find the nearest point xi' = x{rji,rj 2 ) on S to Xg using the 
Newton-Raphson method. (jji,rj 2 ) are the local orthonormal coordinate system on the 
element S. 

• 2. Determine the relative position of the source projection in the local coordinate 
system. 

• 3. Calculate the unit normal ng to the curved element 5 at x^ 

• 4. Approximately project the curved element S on the polygon S in the plane tangent 
to S' at xi' 

• 5. Determine the geometry of the projected quadrilateral S and four triangular 
regions defined by joining x^ to the four corners of S, and then determine the linear 
mapping matrices Lj. 

• 6. Introduce polar coordinates {p, 9) in S, centered at Eq. (|3.2I) becomes: 


I = 





(3.3) 
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where J is the Jacobian of the mapping from Cartesian coordinates on S to curvilinear 
coordinates on S. 

• 7. Apply radial variable transformation: R{p) in order to weaken the near-singularity 
due to which is essentially related to the radial variable p only: 

• 8. Apply the angular variable transformation (0), with hj being for each triangle j 
the distance between and the foot of the perpendicular from 5c7 and the edge of S and 
Uj the angle at of the triangle j: 


t{9) 


hi loff ( l + sin(6>-aj) \ 
2 \1 — sin (0 — Uj)) 


(3.4) 


in order to weaken the angular near-singularity which arises from Pmax{9) when Xg is 
near the edge of the polygon S. This uses the fact that: 

— = ^ - ^j) ^3 5) 

dt Pmax{d^ hj 


• 9. Use Gauss-Legendre’s formula to perform numerical integration of (13.31) in the 
transformed variables R and t : 

. dR (3.6) 


I = 


j 

Jt( 


dt 


rR{Prr.aRe)) J-J^ 


t(0) PmaxiP) 7/{(o) 


dR 


3.3. Determination of physical quantities 

Local surface elevation and its corresponding depth-integrated potential energy are ob¬ 
tained from the free surface data, extracted from the data using the same 3rd order spline 
interpolation polynomials as used in the BEM simulation ('section 12.21) . Velocities and 
the associated depth-integrated kinetic energy are calculated using free surface data and 
the interior helds as seen above (section 13.21) . 

The crest and trough locations are tracked using a two-stage detection algorithm: 

First, extrema and the zero-crossing positions are calculated directly from the dis- 
cretising spline polynomials. These positions are then linked together between time steps 
to form characteristics of the motion. All physical quantities are tracked along these 
characteristic lines. 


4. Results and Discussion 

The numerical simulations undertaken are summarised in Table [T] 







Identifier 

Class 

Number 

of 

Waves 

Approx. 

points 
per wave 

Paddle 

amplitude 

A 

X.focUS 


Cpo 

C3N5A0.05 

3 

5 

8 

0.05 

2D 

1.1014 

1.2878 

C3N5A0.3 

3 

5 

8 

0.3 

2D 

1.1171 

1.22878 

C3N5A0.508 

3 

5 

16 

0.508 

2D 

1.2382 

1.2124 

C3N5A0.511 

3 

5 

16 

0.511 

2D 

1.2499 

1.2028 

C3N5A0.513 

3 

5 

8 

0.513 

2D 

1.2024 

1.2156 

C3N5A0.514 

3 

5 

8/16 

0.514 

2D 

1.2416 

1.1978 

C3N5A0.516 

3 

5 

8/16 

0.516 

2D 

1.2380 

1.2022 

C3N5A0.518 

3 

5 

8/16 

0.518 

2D 

1.2126 

1.2023 

C3N5A0.519 

3 

5 

8/16 

0.519 

2D 

1.2275 

1.1882 

C3N5A0.53 

3 

5 

8 

0.53 

2D 



C3N5A0.56 

3 

5 

8 

0.56 

2D 



C3N5A0.32X10 

3 

5 

8 

0.32 

10 

1.3940 

1.1901 

C3N5A0.33X10 

3 

5 

8 

0.33 

10 

1.2538 

1.2290 

C3N5A0.34X10 

3 

5 

8 

0.34 

10 

1.2472 

1.2250 

C3N5A0.35X10 

3 

5 

8 

0.35 

10 

1.2568 

1.2202 

C3N5A0.36X10 

3 

5 

8 

0.36 

10 

1.2735 

1.2168 

C3N7A0.41 

3 

7 

8 

0.41 

2D 

1.2925 

1.3141 

C3N7A0.42 

3 

7 

8 

0.42 

2D 

1.2974 

1.3155 

C3N7A0.43 

3 

7 

8 

0.43 

2D 

1.2964 

1.3233 

C3N7A0.44 

3 

7 

8 

0.44 

2D 

1.2981 

1.3592 

C3N7A0.45 

3 

7 

8 

0.45 

2D 

1.3075 

1.3082 

C3N7A0.46 

3 

7 

8 

0.46 

2D 

1.3029 

1.3102 

C3N7A0.47 

3 

7 

8 

0.47 

2D 

1.3002 

1.3505 

C3N7A0.48 

3 

7 

8 

0.48 

2D 

1.3156 

1.3085 

C3N7A0.49 

3 

7 

8 

0.49 

2D 

1.3430 

1.3290 

C3N7A0.5 

3 

7 

8 

0.5 

2D 

1.2818 

1.2593 

C3N9A0.42 

3 

9 

8 

0.42 

2D 

1.4651 

1.4521 


o 

Comments 


Small steepness case 
Medium steepness case 

Marginal high resolution recurrence case 

Marginal high resolution breaking case 

Marginal medium resolution recurrence case 
Marginal medium resolution breaking case 
Breaking onset too close to paddle to obtain 
reliable values for Co and Cp^ 

Breaking onset too close to paddle to obtain 
reliable values for Co and Cp^ 

Marginal recurrence case 
Marginal breaking case 


Marginal recurrence case 
Marginal breaking case 
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C3N9A0.43 

3 

9 

8 

0.43 

2D 

1.4771 

1.4942 


C3N9A0.44 

3 

9 

8 

0.44 

2D 

1.4836 

1.5000 


C3N9A0.45 

3 

9 

8 

0.45 

2D 

1.4904 

1.5054 


C3N9A0.46 

3 

9 

8 

0.46 

2D 

1.4938 

1.4551 


C3N9A0.47 

3 

9 

8 

0.47 

2D 

1.5229 

1.5265 

Marginal recurrence case 

C3N9A0.48 

3 

9 

8 

0.48 

2D 

1.6005 

1.4276 

Marginal breaking case 

C3N9A0.49 

3 

9 

8 

0.49 

2D 

1.2659 

1.4044 



Table 1: Summary of computational experiments and their refer¬ 
ence velocities cq and Cp„. 


method 

local 
average 
dispersion plot 


Wo Awo Interpolated wq ko Afco Interpolated ko cq = ^ Interpolated co 

5.3833 0.4141 5.2393 4.1493 0.5187 4.2357 1.2974 0.2020 11405 

5.3833 0.4141 5.2238 4.1490 0.519 4.1793 1.2974 0.2020 1.2499 

5.3833 0.4141 4.1493 0.5187 1.2974 0.2020 

Table 2: Determination of the reference value of the linear veloc¬ 
ity Cq. This table shows the convergence of these values cq using 
3 different methods, local, on averaged data or using the whole 
dispersion plot, in the marginal recurrence C3N5A0.511 case. Aujq 
and Afco are the spectral resolutions around the maximum. 
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Figure 4. C3N5 wave profile with the definitions of the classical and local quantities as well 

as associated polar coordinates. 


4.1. Crest and trough slowdown 

4.1.1. Local quantities 

The local behaviour of the wave kinematics calls for precise definitions of wave geometry 
(FigurelH: amplitude a, height H, wavelength A and crest half wavelength Ac. The local 
crest steepness Sc is defined as Sc = a-^. An analogous definition exists for the local 
trough steepness St- 

Figure [5] shows the relationship between the Stokes steepnes s ak = a nd the local 
crest steepness Sc calculated using a 5*^order Stokes wavetrain (|Fentonlll985l) . It is noted 
that the classical Stokes limit of ak = 0.42 becomes a limiting steepness of Sc = 0.71 
when expressed in terms of the local crest steepness. 


4.1.2. Determination of the reference wave speed cq 

For robust comparison, the reference wave speed cq is required and needs to be precisely 
and carefully determined. A spectral analysis of both time and space wave signals is 
undertaken to be able to estimate the true reference wave quantities rather than those 
based on the paddle signal. 

Figure [S] shows the amplitude spectrum of a representative paddle signal. The char¬ 
acteristic generating frequency ujp is systematically different from the power spectrum 
peak value ujq (vertical black line) obtained by a FFT analysis. Consequently, each in¬ 
vestigated case has to be analysed in space and time to determine the reference power 
spectrum peak wavenumber kg and peak frequency ojq. 

Several methods were used to determine the corresponding values of cq = ^ and 
Cpo = and these show convergence towards identical values for the same simulation 
case. The linear velocity cq has been estimated without making any hypothesis about 
the (deep-water) dispersion relationship, using the peak of the frequency spectrum wq 
and the peak of the wave number spectrum ko, giving the linear celerity from the ratio 
Cq. The deep water dispersion relationship can be used to determine Cp^, using the peak 
frequency wq- 

Table[2]shows an example of three different methods to obtain cq using a local method, 
an average method and a global method for the 2D marginal breaking C3N5A0.511 case. 

The local method value has been found by computing the spatial FFT at a random 
time and a FFT time series centered on the position of the peak of the wave group at 
this given time. The average method is based on determining ujq from the ensemble mean 
FFT of the time series produced by 15 virtual height probes spaced equally, between two 
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Figure 5. Relationship between local steepness Sc and Stokes steepness ak numerically 
computed for a 5*" order Stokes wave. 



Figure 6. Displacement spectrum for the Class 3 N5 paddle stroke signal. The damping com¬ 
ponent modulates the signal and the spectral peak value ojo (solid vertical line) differs from 
the generating wavenumber ujp (dashed vertical line). The spectral peak was located accurately 
using a local spline (solid line). 


crest (or trough) maximum events. The peak wavenumber kg was determined from the 
ensemble mean FFT of 512 surface profiles between two crest (or trough) maxima. The 
global method is a 2D FFT of a space-time domain of the simulation to obtain the 
dispersion graph (figure [7]) . 

The values of cq found in Table [2] show a very small sensitivity to the extent of the 
time and space domain used to determine the reference velocity: cq remains almost con¬ 
stant. These frequency and wavenumber estimates are the raw values found by taking 
the maximum values of the spectra at a given resolution. The length of the time series 
and the size of the domain do not provide sufficiently high resolution in the spectral 
domain (indicated by Awq and Afco)j and a refined method was necessary. An interpo¬ 
lating polynomial is constructed using the neighbouring points around the power spectral 
maximum and the new refined maximum is found analytically using the polynomial co¬ 
efficients. These new quantities give converged values of the reference velocities cq and 
Cpo that have been tabulated for all the cases in the table [T] 

Some additional remarks follow on the value of the reference velocity Cp^ based on the 
deep water dispersion relation. Differences between cq and Cp^ have been found, which 
can be explained by the nature of the wave group. Each individual wave in the group 
has a wavelength short enough to be considered in deep water in the simulation, but 
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uj/u)p 


Figure 7. Dispersion plot for the C3N5A0.511 marginal recurrence case. Colours represent the 
spectral energy density and the axes have been normalised by the paddle frequency LOpaddie. and 
the associated linear wavenumber kpaddie = ojpaddie/g- The coordinates of the peak in the graph 
are used to determine the reference linear velocity co = The solid black line represent the 
Stokes linear deep water dispersion relationship 


the wa yelength of the wave-packet e nvelo pe is long enough not to be considered in deep 
water. Longuet-Higgins fc Steward (1964) called this a transitional wave group. These 
differences are secondary and do not affect the results. 

In addition to the information about the spectral peak waves, the dispersion plot shows 
the degree of nonlinearity obtained with the multiple higher frequency “ripples”. These are 
the bound waves of the group, i.e. higher harmonics of the basic free waves constituting 
the wave group. In general, the more harmonics that are present in the dispersion graph, 
the more nonlinear is the whole phenomenon. 

4.2. Space-time characteristics 

This study has been carried out using chirp ed Class 3 wave pack ets as defined above, but 
broader experimental and field studies (see Banner et all ( 2014h l confirm its generality. 

Figure [5] represents the characteristics plot of crests and troughs along the centre slice 
of the 3D numerical wave tank of a representative Class 3 wave packet. Crests, troughs 
and zero-crossings have been detected and followed in this coordinate system, the abscissa 
is the non-dimensional position and the ordinate is the non-dimensional time. In this plot, 
the theoretical linear trajectory is represented by a long-dash line as a reference, crest 
paths are the full lines and trough paths are represented by dashed lines. Trajectories 
in space and time of each crest and trough are clearly seen as nonlinear, with local 
variations. 

This plot exhibits systematic patterns for the crest and trough trajectories during the 
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Crests 

■ ■ Troughs 



Figure 8. Space-time plot of crests (full lines) and troughs (dashed lines) of the C3N5 marginal 
recurrence case. The long-dashed line represents the reference linear velocity. Slowdown region 
occurs where the lines turn anticlockwise. In this plot, a systematic slowdown and reacceleration 
of each crest and trough is visible. 


life cycle of each wave observed. The local slope of these trajectories represents the veloc¬ 
ity of the crest or trough. A trajectory path turning anticlockwise implies slowdown and 
a clockwise rotation represents an acceleration. Each of the paths has a clear systematic 
variation in the slope and slowdown can be observed before a re-acceleration, between 
x/Xp = 1 and x/Xp = 5. These slowdowns occur close to the crest and trough maxima 
of each wave inside the wave packet. 


4.3. Slowdown of each wave, a generic feature. 


Slowd o wn phenomena have been obs e rved and mentioned prev iously (jJohannessen &: Swan 
( 2001 1. .Tohannessen fc Swan (2003), Katsardi fc Swaiil ( 201lll ). but have not been stud¬ 
ied as a phase-resolved wave effect. The following plots address this gap by presenting 
the evolution of the crest (or trough) velocity as a function of the local steepness Sc- 
The curves show the trajectory of the non-dimensional horizontal velocity of each carrier 
wave crest in the wavegroup, normalised by the theoretical linear wavespeed cq defined 
above. Velocities are calculated using the first derivative of a 3rd-order spline fitting the 
space-time trajectory. Results below are presented for crests only for illustration, and 
troughs are not shown because they behave in exactly the same way. 

Figure [5] shows the evolution of the horizontal crest velocity of every crest of a Class 
3, N = 5 marginal recurrent wave packet. The abscissa measures the non-dimensional 
velocity ratio ^ and the ordinate shows the local steepness Sc- Each individual crest 
life-cycle in the group is represented by a line trajectory. A black dot indicates where 
each trajectory starts. 

Two important finding are represented in this plot showing the instantaneous crest 
velocity against the local steepness. The first conclusion, which is one of the new findings 
of this study, is that contrary to inferences based on Stokes wavetrain theory, the crest 
velocity is not strongly dependent on the steepness. Minimum crest velocity in the en¬ 
semble of waves shown in figure [9] approaches the same value systematically, with a local 
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Figure 9. Local steepness plotted against normalised velocity for a C3N5A0.511 marginal 
recurrent wavegroup. The dashed line represents 5th order Stokes celerity for reference. Each of 
the lines represents an individual crest within the group which undertakes an hysteresis cycle, 
showing a slowdown as the local steepness increases. The starting point is indicated by a black 
dot. Spikes seen in the velocity trajectories are artefacts from the tracking algorithm and are 
explained below. 


steepness varying from Sc = 0.33 to Sc = 0.52. Previous knowledge of velocity variation 
of crests is represented by the dashed line, which corresponds to 5*^order Stokes crest 
celerity. 

The second new result in this plot concerns the minimum crest velocity. In this marginal 
recurrent case, the velocities of each crest reduce close to a minimum of 0.8 of the linear 
reference velocity cq associated with this packet. The evolution of velocity is not linear 
with the steepness, and loops can be seen in the curves describing the history of the crest 
velocity, which shows the existence of hysteresis behaviour. It is noted that the minimum 
speed of this ensemble of crests seems to rise a few percent (from ^ = 0.75 to 0.78 in 
this case) when the local steepness ri ses and the strong nonlinear i ties enter. _ 

The s e findings are corroborate d by .lohannessen fc SwaiJ ( 200l[ l , .lohannessen fc Swan 
I 2003h . iKatsardi fc Swan ( 2011 1. These authors show the existence of a slowdown of 
10% in 2D waves and of 7% in the case of 3D waves, with scatter. These slowdowns are 
quantified in respect to a reference linear velocity, called ui , which is not clearly defined. If 
a guess of the reference velocity estimation based on the peak of the frequency spectra is 
made, then the normalising velocity becomes Cpg defined above. With this hypothesis, the 
minimum velocity for the dominant crest in this graph goes from — = 0.78 to — =0.80, 

Cq CpQ 


resulting in a 20% reduction in wave speed. A key question that arises from this analysis 
is whether this effect is local. The above papers argue that the cause of the slowdown 
can be associated with the appearance in the frequency wave spectrum of a peak around 
1.2 to 1.5 of the main wave-packet peak, producing slower waves in accordance with the 
theoretical linear free-wave dispersion relationship. However, the present study offers a 
complementary argument for the cause of the slowdown in section 14.41 below. 

The next figure addresses whether the slowdown is only a local crest and trough effect. 
Figure [TO] shows the velocity of the mid-point between the two zero-crossings which 
adjoin (or span) each crest (or trough) studied. Each zero-crossing pair has the steepness 
attribute of the crest they span. The abscissa is the non-dimensional velocity speed with 
respect to the reference linear phase velocity cq . 

The crest velocity at each crest maximum is represented by a dot, and the bars repre- 
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Figure 10. Zero-crossing mid-point velocity as a function of the max;imum local steepness Sc 
for each crest for a range of 3D C3N5, 2D C3N5 and C3N7 cases, plotted against the crest 
velocity. The velocity at each crest maximum is represented by a dot, and the bars represent 
the extent of the variation of the velocity during each individual crest life-cycle. 


sent the extent of the variation of the crest velocity during each individual crest life-cycle. 
The plot shows that the zero-crossing speed changes dynamically during the wave evo¬ 
lution, and can undertake large relative speed variations of the order of roughly 20%. 
Nonetheless, the speed values when the crest reaches its maximum (dots) are close to the 
linear velocity, within the range O.Ocq — I.Ocq. Consequently, the slowdown is a localised 
crest effect, and does not affect the whole packet uniformly. 

4.4. Leaning, cause of the slowdown. 

This section will show that the cause of the slowdown is geometric. We define the leaning 
as the left/right asymmetry : 


L = 


-1 + 2 


X-w 

Xr-Xi 


(4.1) 


where X represents the crest (trough) position, Xi is the position of the downward 
zero-crossing and Xr is the position of the upward zero-crossing. The leaning L quantifies 
the asymmetry of the geometry of the triangle formed by the two zero-crossings spanning 
the crest (or trough). This parameter varies from —1 when the backward-leaning face 
has a vertical asymptote to +1, when the forward-leaning face is vertical. 

Figure m defines polar coordinates associated with a local crest (or trough). 9 is the 
angle at the mid-point (half-way point between the left and right zero-crossings) between 
the zero water level and the crest. R is the radius, the distance between the crest and the 
middle-point. Using polar decomposition of the velocity, the projection on the x axis of 

d6 

the ortho-radial component (the rotation part) is Cqx = —R — sin(^) and the projection 

at 

dR 

on the x-axis of the radial component is Crx = — cos(^). 

dt 
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Figure 11. Ratio of sum of Ce^ (a;-axis projection of radial polar velocity) plus the zero-crossing 
speed plus c (crest velocity at maximum) against maximum steepness Sc for an ensemble of crests 
for 3D C3N5, 2D C3N5 and C3N7 cases. Each dot represent an individual crest. This graph 
shows the crest slowdown is almost entirely included in the radial relative velocity, implying 
that the rotation (leaning) is the cause of the slowdown. 


Figure [TT] shows the sum of the normalised frame of reference velocity 


idX, 


dt 


and the 2 :-projection of the normalised ortho-radial component of the relative velocity 


= —-^■^sin(0) at the crest maximum on the abscissa against the local steepness 


de 


Sc for C3 N5, C3N7, 2D and 3D wave groups. This plot shows a very strong correlation 
between the rotation component and the crest velocity, and proves that the major part 
of the slowdown resides in the rotation (hence the leaning) of crests about the mid-point 
between the zero-crossings. For all the waves in the group, each individual wave tends to 
lean from forward to backward, reducing the crest velocity to below the expected linear 
crest velocity. 

As an example, figure [T^ summarises several key kinematic and geometric properties 
of the dominant crest of a C3N5 wave group. In these graphs, there is a zero-reference 
time (along the abscissa) when the crest reaches its maximum elevation and the local 
steepness is normalised to 1 at this crest maximum. The maximum local steepness for 
this wave is Sc = 0.512. 

This plot represents the entire history of a single carrier wave, with the first panel 
showing the normalised steepness and the second panel showing the leaning parameter 
L. The third panel shows the average zero-crossing speed ^ normalised by the linear 
reference velocity, the fourth panel shows the crest velocity ^ and the last panel is 


proportional to the leaning velocity The graph reproduces the evolution history 

of the dominant crest travelling from the rear of the wave group towards the front of the 
group. 

The wave peaks at tu)p/2'n' = 0. There is a strong correlation between the crest leaning 
backwards, which produces a negative time derivative of L, and the crest slowdown, as 
seen between tujpl2T: = — 1 and ta;p/27r = 1. During the same time interval, the average 
speed of the mid-point of the adjacent zero-crossings is almost constant, with a value 
close to the linear reference velocity. Some peaks appear in the velocity and the time 
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Figure 12. Time evolution of the normalised wave properties tracking the dominant crest of 
a C3N50.511 marginally-recurrent wave group. The first panel is the local crest steepness, the 
second panel is the leaning factor, the third panel is the mid-point zero-crossing velocity, the 
fourth panel is the crest velocity and the bottom panel is the time derivative of the leaning. Note 
how well correlated are the quantities at crest maximum (t=0): Leaning L is 0, meaning the 
wave is symmetric, zero-crossing speed is close to 1.1, c is around 0.8, and the time derivative 
of L indicates the symmetry is changing steadily backwards. 


derivative of L, and are associated with a strong forward jump in L at tujpl2'K = —0.9. 
This is just an artefact of the crest-tracking process and is addressed in the section (14.51) . 

One important hnding is that L = 0 always when the crest (or trough) reaches its 
maximum, implying that the wave crest (or trough) shape is spatially symmetric at 
this instant. This leaning property could potentially be used to detect crest and trough 
maxima in a wave field. These findings imply that the crest slowdown is a strong kinematic 
effect, easily understood via the composition of velocity law: the frame of reference is the 
average position of the two adjacent zero-crossings travelling at roughly the linear phase 
velocity, and the local velocity in this frame of reference is associated with the leaning 
forward-to-backward tendency of the crest, inducing a negative relative crest velocity. 
The composition of the frame velocity and the backward leaning reduce the absolute 
crest velocity to below the originally-expected linear phase velocity. Other fluctuations 
in the absolute velocity confirming this analysis can be seen just after ta;p/27r = 1, where 
a positive velocity of the frame of reference increases the absolute velocity above the 
linear level. 

Figure [13] presents multiple wave profiles stacked in order to form a space-time diagram 
of the typical evolution along the centreline of a C3N5 maximum recurrent wave group. 
The axes have the non-dimensional position as the abscissa and non-dimensional time 
as the ordinate. The wave profiles are colour-coded with the value of the leaning. This 
figure summarises the behaviour described above: both crests and troughs undertake 
a systematic evolution throughout their life cycle. It can be generalised as a generic 
behaviour as follows: a crest (or a trough) begins its life leaning forward (red shading), 
peaks (the geometry of the wave becomes symmetric (green shading)) and decays leaning 
backwards (blue shading). Some abrupt jumps towards the front of the wave (seen as 
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Figure 13. Space-time plot of the evolution of the surface profile towards the C3N5 group 
maximum. Spatial profiles are stacked in time and colour-coded to highlight the leaning dynam¬ 
ics. Green dots indicate crest positions, red dots indicate troughs, bold black stars are crests 
and troughs below the threshold of detection and black dots show zero crossings. Leaning is 
colour-coded between zero-crossings: red (forward-leaning) to blue (backward-leaning). The lin¬ 
ear Stokes velocity co is shown as a blue long-dashed line. Each crest and trough undertakes a 
leaning cycle from leaning forward, peaking at symmetry and then leaning backwards. When the 
crests/troughs slow down (trajectories turn anticlockwise), there is a clear change of symmetry 
when the extrema transition from leaning forward to leaning backwards. 



sudden changes of shading colour) are also present and are explained in the following 
section (I43]| . 

4.5. Apparent scattering and detailed geometric analysis: crestlets or riding waves? 

This section explains the apparent velocity peaks found in the velocity plots (see Fig [12] 
around t = —0.8, also apparent in Fig [HI and Fig [T51) . Figure [T^ is a space-time plot of the 
same wave packet shown in Figure [T51 The abscissa represents the position, the ordinate 
represents the time coordinate and the depth represents the non-dimensional free surface 
elevation of a slice along the centreline of the numerical wave tank. The plot is a time 
sequence of wave profiles. 

In this section we explain the limitation of the crest-tracking algorithm and the appar¬ 
ent velocity spikes. The red arrows in Figure 14 identify where the phenomenon occurs. 
The idea of a crest or a trough being a single maximum or minimum between two zero- 
crossings (which is an arbitrary reference) refers to a Stokes wave. Natural unsteady 
waves sometimes have multiple crests and troughs b etween the same two zero-crossings, 
which evolve dynamically. This effec t can be seen in Johannessen ( 201CJ) in their obser¬ 
vational surface elevation plots. Also, Bateman et all ( 2012 1 describes a numerical model 
results where crestlets and leaning in the plots can be seen in their figures 8 and 9. 
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Figure 14. Space-time plot of the centreline of the 2D C3N5 marginal recurrent case. 2D profiles 
are stacked in the time direction. Spikes found in celerity and leaning values are explained 
by the crest-trough tracking algorithm. A crest/trough has an evolution characterised by a 
succession of ’crestlets’ which appear on the forward face and propagate backwards, and will 
become dominant eventually. Each time the forward crestlet becomes dominant, the detected 
crest position suddenly jumps forward creating these spikes, as highlighted by the red arrows. 


During the growth of a crest, there is an interchange of these small crests or troughs, 
the first one evolving and then the next one appearing from the front of the packet to 
become temporarily the dominant crest, and then the next one appears. These crests are 
initiated, in the relative frame of reference of the wave group, at the front of the group 
and then travel rearwards. During this crest exchange, which is rapid (roughly a tenth 
of a time unit over a tenth of the reference length), the tracking system jumps from the 
decaying crest to the next growing crest, sometimes going through a ’flattened’ area of 
the crest. This creates a discontinuity in the detected crest position, as it is now located 
further towards the front of the group. This produces the apparent high velocity spike. 
It should be noted that between each of these peaks, the celerity of these "crestlets" 
decreases, back to O.Scq. 


The cause of thes e leaning effects is not yet clear, but possible insight is provided by 
Bettini et al\ 1 1983ll who a discuss the radiative tail of a soliton. The radiative tail is the 


train of trailing waves produced by a propagating wave group. In our application, these 
forward jumps of the crests (or troughs) align on the characteristic plot (Figure [5]). The 
average slope gives a velocity one half of the linear velocity, this being the definition of 
the radiative tail in Bettini’s article. These crestlets (or troughlets) seem to propagate 
backwards. This effect is relative because their velocity is one half of the velocity of 
the wave crests (or troughs), which causes the illusion. This energy leakage creates the 
radiative tail. 





































22 


X. Barthelemy, M.L. Banner, W.L. Peirson, F. Dias, M. Allis 



Figure 15. Ratio of computed velocity profiles of a catalogue of cases beneath wave crests to 
the equivalent 5*^ order Stokes velocities, as a function of normalised depth. The computed 
values have been corrected for the Eulerian drift 


5. Subsurface kinematics 

5.1. Velocity profiles 

This section provides further details on the subsurface motion by investigating the veloc¬ 
ity profiles. Figure [15] shows the ratio of the velocity profiles computed from the NWT 
simulation and the 5th order Stokes velocity profiles for different cases. Shown here are 
results below the local crest maxima of 2D C3N5 and C3N9 breaking and non-breaking 
(low steepness to marginal recurrence) , as well as some 3D cases. The Stokes waves used 
for comparison were fitted according to height H and using the neighbouring troughs 
to define the wavelength. Stokes velocity profiles have been corrected for the Eulerian 
current by matching velocities at the bottom of the tank. 

Different zones can be identified. The first zone, around kz = —4, shows a relative 
ratio which can grow up to 1.8. While this value appears substantial, the actual velocity 
differences are not very large because the absolute values of the velocities are small, due 
to the quasi-exponential attenuation of the velocity with depth. 

The second zone, above kz = 0, is more important as it is close to the free surface. Here 
the ratio varies between 1.2 for maximum recurrence cases and 1.7 when the wave breaks. 
This change is significant and highlights an important difference between the predicted 
Stokes orbital velocity and the actual orbital velocity computed by the simulation. The 
consequence can be important when estimating the impact of a wave on a structure 
where the drag is proportional to the square of the velocity. The Stokes fit systematically 
underestimates the particle crest velocity by at least 20% for very steep waves, suggesting 
a near-halving of the impulse loading. 
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Figure 16. Depth-integrated energy partitioning averaged over the wave group for multiple 2D 
C3N7 cases. The evolution of breaking cases is shown by the bottom dot-dashed line ending at 
t = 14.2, and the second (solid) line from the bottom and the third (solid) line from the bottom, 
both terminating at t = 16.3. The ratio of potential energy to total energy of these unsteady, 
nonlinear wave groups shows that linear equipartitioning is almost conserved at this scale. Small 
variations are evident, with breaking groups tending to have slightly more kinetic energy as the 
ratio in these cases decreases to 0.49 

6. Energy partitioning 

One of the main properties commonly assumed for waves is equipartitioning of the 
mean kinetic energy (KE) and potential energy (PE). With this hypothesis, meaurements 
can be made that only measure wave heights (hence PE), from which the kinetic energy 
is deduced. This equipartitioning is theoretically true for steady linear wave trains, but 
does it hold for a nonlinear unsteady wave group? This section investigates the energy 
partitioning in a global sense, first analysing the whole wave group and then studying 
locally what occurs at the scale of the individual crest and troughs. 

6.1. Global partitioning 

Having proven in the simulation validation section that energy conservation was re¬ 
spected, further investigation was undertaken concerning the energy partitioning for the 
whole wave group. Eigure (fTHl) shows the depth-integrated energy partition for the 
whole wave group with respect to time, for the C3N7 breaking and non-breaking cases. 
Periods when the whole group was not in the integration domain have been blanked out. 
This figure shows that global equipartitioning holds with a small excess in the kinetic 
energy at the crest maximum. The ratio oscillates globally close to = 0.5, showing 
that the total energy (TE) is almost equally divided between potential energy (PE) and 
kinetic energy (KE). Breaking cases have the maximum deviation from equipartitioning 
with ~ 49% of the energy residing in the PE. The energy flux necessary to produce the 
jet formed as the breaking event initiates has its source in the kinetic energy. 

6.2. Local energy partitioning at crests and troughs 
The above section showed how the depth-integrated energy partitioning still holds for 
the whole unsteady wave group. We now examine the local behaviour of the depth- 
integrated energy partitioning at both crest and trough maxima. Eigure flTl shows the 
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Figure 17. Depth-integrated energy partitioning for the dominant crests (C3N5 are circles, 
C3N7 are stars, C3N9 are pluses, C3N10 are x-crosses and 3D C3N5 are pentagons) and troughs 
(C3N5 are inverted triangles, C3N7 are triangles, C3N9 are left-pointing triangles, C3N10 are 
right-pointing triangles and 3D C3N5 are diamonds) at local group maximum for a catalogue 
of cases. The dashed lines represent this ratio for a 5th order Stokes wave. A clear deviation 
from linear theory energy partitioning is seen, with 2D troughs clustering around a ratio of local 
depth-integrated potential energy to local depth-integrated total energy of 0.7, and 3D troughs 
around 0.71. Also, 2D and 3D crests cluster from 0.6 to 0.64. This ratio drops below 0.6 at the 
initiation of breaking crests, signalling an excess of KE. Small amplitude waves have the same 
ratio of 0.67 for crests and troughs. 


energy partition plotted against local steepness Sc for 2D and 3D breaking for the 
maximum recurrent wave packet cases C3N5, C3N7, C3N9 and C3N10. The dashed lines 
represent this ratio for a 5th order Stokes wave. Triangles represent the values for the 
troughs and circles show the values for the crests. Remarkably, our results show that 
equipartitioning breaks down locally at crest and trough maxima and the ratio stays 
almost constant, 0.6 to 0.62 for the crests and 0.7 for the troughs. Locally, the Stokes 
wave does not conform to the equipartitioning, but behaves differently from the unsteady 
group waves. The ratio for the unsteady wave groups matches the Stokes wave for 
very small steepness and is close to 0.67, but diverges from this level at higher steepness. 
The notable difference between the marginal recurrent case and the breaking case is the 
breaking crest acquires proportionally more kinetic energy with the partition ratio falling 
below 0.6. It is again evident that the source of the breaking jet resides in the kinetic 
energy, while the wave height (hence the potential energy) remains almost constant. 


7. Conclusions and recommendations 

This computational study of highly nonlinear gravity water waves unsteady wave 
groups has revealed new insights into the behaviour of kinematic and geometric prop¬ 
erties, as well as the associated energetics. In particular, our findings on generic crest 
slowdown obtained from these numerical chirped wave group simulations have been con- 
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firmed observationally for chirped as well as for bimodal wave groups in our laboratory 
wave basin experi ments, as well as for deep water wave groups measured in an ocean 


tower experiment ( Banner et al. ( 2014ll l. 


When the theoretical constraint of steadiness is relaxed, additional asymmetric degrees 
of freedom are observed in the wave shape. In particular, the kinematics reveal a forward- 
to-backward leaning cycle for each individual crest and trough inside the wave group, 
resulting in a local crest or trough slowdown of around 20% of the reference linear velocity 
estimated from a dispersion plot. This is a geometric effect, as the velocity of the zero- 
crossings remains close to the linear phase velocity and the observed relative velocity 
variations can be explained by the geometric changes in the waveform. Velocities profiles 
under the crest show a large difference between the Stokes-fit velocity profile and the 
computed velocity profile. Computed orbital velocities under the crest can be signihcantly 
higher than for the Stokes wave of corresponding height, and the ratio between them 
varies between 1.2 and 1.8 for the steepest cases. These values should be taken in account 
when designing structures close to the surface of the ocean where the Stokes values are 
often assumed. The associated energetics have been carefully monitored and support, 
within 1%, equipartitioning of the depth-integrated energy on the scale of the wave 
group, although at breaking onset, a slight excess of kinetic energy occurs. This depth- 
integrated energy equipartitioning breaks down when investigated locally at crest and 
trough maxima during the evolution of the wave group. These show an energy partitioning 
of around 69% for the potential energy at trough maxima and a partitioning range from 
60% to 64% of the potential energy at crest maxima. The crest energy partitioning 
at breaking drops below 60% for the potential energy, indicating the initiation of the 
breaking jet. 

Overall, these findings contribute new insights into how actual ocean wave groups 
propagate when the academic constraint of steadiness is relaxed. Aside from its intrinsic 
interest, the generic crest slowdown explains the reduced speed of initiation of breaker 
fronts that has been reported occasionally in the literature and is central to assimilating 
whitecap data accurately into spectral sea-state forecast models. Also, parameterisations 
of air-sea fluxes of momentum and energy, which depend on the square and cube of the 
sea-surface fluid velocity, may also be modified appreciably, as may the impact loading 
on offshore structures. 
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